Part 2 Fire Selection

2.1 Set Up

2.1.1 Libraries

library(tidyverse)
library(terra)
library(sf)
library(mapview)
library(raster)
library(rgeos)
library(lubridate)
library(ggplot2)
library(exactextractr)
library(gridExtra)
library(knitr)
library(rasterVis)
library(RColorBrewer)
library(spData)
library(forcats)
library(cowplot)
library(rgeos)

2.1.2 USDA National Forest Type Group Dataset

Conifer Forest Type Groups: Douglas-Fir, Fir-Spruce-Mountain Hemlock, Lodgepole Pine

# forest type groups and key
conus_forestgroup <- raster('data/forest_type/conus_forestgroup.tif')
forest_codes <- read_csv('data/forest_type/forestgroupcodes.csv')

# set crs
crs = crs(conus_forestgroup)

2.1.3 EPA level-3 Ecoregions

Canadian Rockies, Idaho Batholith, Middle Rockies, Columbian Mountains - Northern Rockies

# level 3 ecoregions
l3eco <- st_read('data/ecoregion/us_eco_l3.shp') %>%
  st_transform(., crs=crs)

# select northern rocky mountains from level3 ecoregions
eco_select <- l3eco %>% 
  filter(NA_L3NAME %in% c('Canadian Rockies','Columbia Mountains/Northern Rockies','Middle Rockies','Idaho Batholith'))

2.1.4 Mapping

2.1.4.1 Ecoregions

# mapview
palette <- brewer.pal(18,"YlGn")
palette[1] <- rgb(255, 255, 255, maxColorValue=255, alpha=1)
mapview(eco_select,na.color=palette[1],zcol = "US_L3NAME",layer.name = "Level-3 Ecoregion", legend=TRUE)

2.1.4.2 Forest Type Groups

# convert raster values to factors
forestgroup_eco <- crop(conus_forestgroup,eco_select) %>% 
  mask(.,eco_select) 

forestgroup_eco[forestgroup_eco %in% c(120,180,240,300,320, 360, 370, 400, 500,700,900,920,950)] <- 1

# add a labels for forest type code 
group_levels <- levels(forestgroup_eco %>% as.factor())[[1]]
group_levels[["Forest Type"]] <- c("Unforested","Other","200: Douglas-fir","220: Ponderosa Pine","260: Fir/Spruce/Mountain Hemlock","280: Lodgepole Pine")
# group_levels[["forest_type"]] <- c("Unforested","120: Spruce/Fir","180: Pinyon/Juniper","200: Douglas-fir","220: Ponderosa Pine","240: Western White Pine","260: Fir/Spruce/Mountain Hemlock","280: Lodgepole Pine","300: Hemlock/Sitka Spruce","320: Western Larch","360: Other Western Softwood","370: California Mixed Conifer","400: Oak/Pine","500: Oak/Hickory","700: Elm/Ash/Cottonwood","900: Aspen/Birch","920: Western Oak")

levels(forestgroup_eco) <- group_levels

# mapview
mapview(forestgroup_eco, col.regions=palette,na.color=palette[1],legend=TRUE,layer.name = "Forest Type")

2.2 Define Fire Parameters

2.2.2 Group Adjacent Fires

# function to group adjoining fire polygons to ensure contiguous high-severity patches across MTBS events
group_fires <- function(mtbs_year) {

  # join the polygons with themselves, and remove those that do not join with any besides themselves
  combined<- st_join(mtbs_year, mtbs_year, join=st_is_within_distance, dist = 180, left = TRUE,remove_self = TRUE) %>% 
    drop_na(Event_ID.y)%>% 
    dplyr::select(Event_ID.x,Event_ID.y)
  
  if(nrow(combined)>=1){ # if there are overlaps for this years fires...
    
    # partition data into that that has overlap, and that that does not
    overlap <- mtbs_year %>%
      filter(Event_ID %in% combined$Event_ID.x)
    no_overlap <- mtbs_year %>%
      filter(!(Event_ID %in% combined$Event_ID.x))
    
    print(paste0("there are ",nrow(overlap)," overlapping polygons"))
    
    # join all overlapping features, and buffer to ensure proper grouping
    overlap_union <- st_union(overlap) %>%
      st_buffer(190)
    
    # break apart the joined polygons into their individual groups
    groups <- st_as_sf(st_cast(overlap_union ,to='POLYGON',group_or_split=TRUE)) %>%
      mutate(year = mean(mtbs_year$year),
             Fire_ID = str_c("Fire_",c(1:nrow(.)),"_",year)) %>%
      rename(geometry = x)
    
    print(paste0("polygons formed into ",nrow(groups)," groups"))
    
    # join back with original dataset to return to unbuffered geometry
    grouped_overlap <- st_join(overlap,groups,left=TRUE)
    
    # arrange by the new grouping
    joined_overlap_groups <- grouped_overlap %>%
      group_by(Fire_ID) %>%
      tally()%>%
      st_buffer(1) %>%
      dplyr::select(Fire_ID) %>%
      mutate(year = mean(mtbs_year$year))
    
    # add new ID to the freestanding polygons
    no_overlap_groups <- no_overlap %>%
      mutate(Fire_ID = str_c("Fire_",nrow(groups)+c(1:nrow(no_overlap)),"_",year)) %>%
      dplyr::select(Fire_ID,year)
    
    # join the new grouped overlap and the polygons without overlap
    fires_export <- rbind(joined_overlap_groups,no_overlap_groups)
    return(fires_export)
    
    } else { # if there are no overlaps for this year...
      
      print("no overlapping polygons")
      
      fires_export <- mtbs_year %>%
        mutate(Fire_ID = str_c("Fire_",c(1:nrow(.)),"_",year)) %>%
        dplyr::select(Fire_ID,year)
      
      return(fires_export)
  }
}
# group adjacent polygons within each fire year
fires_88 <- group_fires(mtbs_select %>%  filter(year == 1988))
## [1] "there are 22 overlapping polygons"
## [1] "polygons formed into 7 groups"
fires_89 <- group_fires(mtbs_select %>%  filter(year == 1989))
## [1] "there are 2 overlapping polygons"
## [1] "polygons formed into 1 groups"
fires_90 <- group_fires(mtbs_select %>%  filter(year == 1990))
## [1] "there are 2 overlapping polygons"
## [1] "polygons formed into 1 groups"
fires_91 <- group_fires(mtbs_select %>%  filter(year == 1991))
## [1] "no overlapping polygons"
# join each fire year
mtbs_grouped <- rbind(fires_88,fires_89,fires_90,fires_91)%>%
  mutate(area_ha = as.numeric(st_area(geometry))/10000,
         area_acres = area_ha*2.471)

2.2.3 Select Fires by Ecoregion and Forest Type

# assign ecoregion and proportions of forest type to each fire polygon
fires_join <- st_join(mtbs_grouped,eco_select,join=st_intersects,left=FALSE,largest=TRUE) %>% 
  left_join(., exact_extract(conus_forestgroup,mtbs_grouped, append_cols = TRUE, max_cells_in_memory = 3e+08, 
                             fun = function(value, coverage_fraction) {
                               data.frame(value = value,
                                          frac = coverage_fraction / sum(coverage_fraction)) %>%
                                 group_by(value) %>%
                                 summarize(freq = sum(frac), .groups = 'drop') %>%
                                 pivot_wider(names_from = 'value',
                                             names_prefix = 'freq_',
                                             values_from = 'freq')}) %>%
              mutate(across(starts_with('freq'), replace_na, 0)))
 
# remove unnecessary columns, cleanup names
# filter to ensure fire polygons are at least 25% type of interest
fires <- fires_join %>% 
  dplyr::select("Fire_ID","year","area_ha","area_acres","US_L3NAME","freq_0","freq_200","freq_220","freq_260","freq_280") %>% 
  rename("ecoregion" = "US_L3NAME",
         "freq_df"="freq_200",
         "freq_pp"="freq_220",
         "freq_fs"="freq_260",
         "freq_lpp"="freq_280") %>% 
  mutate(freq_allother = 1-(freq_0 + freq_df+freq_pp+freq_fs+freq_lpp),
         freq_forested = 1- freq_0,
         freq_ideal = freq_df+freq_fs+freq_lpp)%>% 
  mutate(across(starts_with('freq'), round,2))%>% 
  filter(freq_ideal > 0.25)

2.2.4 Select Fires by Burn Severity

# import all mtbs rasters via a list
rastlist <- list.files(path = "data/mtbs", pattern='.tif', all.files=TRUE, full.names=TRUE)
allrasters <- lapply(rastlist, raster)
names(allrasters) <- str_c("y", str_sub(rastlist,22,25))

# create empty dataframe
severity_list <- list()

# loop through mtbs mosasics for 1988-1991
# extract mtbs burn severity raster for all selected fires
# calculate burn severity percentages for each fire
for (i in names(allrasters)){
  mtbs_year <- allrasters[[i]]
  fire_year <- filter(fires, year==str_sub(i,2,5)) 
  raster_extract <- exact_extract(mtbs_year,fire_year, max_cells_in_memory = 3e+09,coverage_area=TRUE)
  names(raster_extract) <- fire_year$Fire_ID 
  
  output_select <- bind_rows(raster_extract, .id = "Fire_ID")%>%
    group_by(Fire_ID , value) %>%
    summarize(total_area = sum(coverage_area)) %>%
    group_by(Fire_ID) %>%
    mutate(proportion = total_area/sum(total_area))%>% 
    dplyr::select("Fire_ID","value","proportion") %>% 
    spread(.,key="value",value = "proportion")
  
  severity_list[[i]] <- output_select
}

# combine extracted raster datasets
severity_df <- do.call(rbind, severity_list) 

# join burn severity % to fires polygons
# filter dataset for 500 acres high severity
fires_severity <- left_join(fires,severity_df,by="Fire_ID")%>% 
  rename(noburn= "1",lowsev = "2", medsev = "3", highsev = "4",regrowth = "5", error = "6") %>% 
  dplyr::select(- "NaN",-"regrowth",-"error") %>% 
  mutate(highsev_acres = area_acres*highsev)%>% 
  filter(highsev_acres > 500)

2.2.5 Clean Up Dataset

# get the most common forest type within each polygon
fires_select <- fires_severity %>%
  left_join(.,exact_extract(conus_forestgroup,fires_severity, 'mode', append_cols = TRUE, max_cells_in_memory = 3e+08)) 

fires_select$mode <- as.factor(fires_select$mode)

fires_select <- fires_select %>% 
    mutate(fire_foresttype = case_when(mode==200 ~ "Douglas-Fir",
                                       mode==220 ~ "Ponderosa",
                                       mode==260 ~ "Fir-Spruce",
                                       mode==280 ~ "Lodegepole Pine",
                                       TRUE ~ "Other"),
           Fire_ID = str_c("Fire_",c(1:nrow(.)),"_",year))
# join the grouped fires back to original mtbs boundaries
fires_mtbs <- st_join(mtbs_select,fires_select,left=FALSE,largest=TRUE) %>% 
  filter(year.x==year.y)%>% 
  dplyr::select("Event_ID","Incid_Name","Fire_ID","Ig_Date","year.y","state","BurnBndAc","ecoregion") %>% 
  rename(year= year.y)

2.3 Final Fire Dataset

2.3.1 Dataset Overview

full_dataset <- fires_select %>% 
  st_drop_geometry() %>% 
  dplyr::select(Fire_ID,year,ecoregion,fire_foresttype,area_acres,highsev) %>% 
  mutate(highsev = round(highsev,2),
         area_acres = round(area_acres,0))  

kable(full_dataset, 
      align = 'c',
      padding = 1,
      col.names = c("Fire ID", "Year", "Ecoregion", "Majority Forest Type","Area (acres)", "High Severity %"),
      caption = "High-Severity Conifer-Dominated Fires 1988-1991")
Table 2.1: High-Severity Conifer-Dominated Fires 1988-1991
Fire ID Year Ecoregion Majority Forest Type Area (acres) High Severity %
Fire_1_1988 1988 Middle Rockies Fir-Spruce 342005 0.27
Fire_2_1988 1988 Middle Rockies Lodegepole Pine 777690 0.28
Fire_3_1988 1988 Middle Rockies Lodegepole Pine 448911 0.21
Fire_4_1988 1988 Idaho Batholith Douglas-Fir 5651 0.16
Fire_5_1988 1988 Idaho Batholith Fir-Spruce 11945 0.23
Fire_6_1988 1988 Idaho Batholith Douglas-Fir 50666 0.07
Fire_7_1988 1988 Middle Rockies Lodegepole Pine 167870 0.50
Fire_8_1988 1988 Idaho Batholith Douglas-Fir 25889 0.02
Fire_9_1988 1988 Idaho Batholith Douglas-Fir 8312 0.17
Fire_10_1988 1988 Canadian Rockies Lodegepole Pine 42492 0.52
Fire_11_1988 1988 Idaho Batholith Fir-Spruce 45075 0.43
Fire_12_1988 1988 Middle Rockies Lodegepole Pine 5633 0.24
Fire_13_1988 1988 Idaho Batholith Douglas-Fir 4962 0.25
Fire_14_1988 1988 Middle Rockies Other 35864 0.46
Fire_15_1988 1988 Idaho Batholith Fir-Spruce 19499 0.29
Fire_16_1988 1988 Idaho Batholith Fir-Spruce 5626 0.09
Fire_17_1988 1988 Idaho Batholith Fir-Spruce 12746 0.40
Fire_18_1988 1988 Middle Rockies Other 13108 0.48
Fire_19_1988 1988 Idaho Batholith Fir-Spruce 7241 0.27
Fire_20_1988 1988 Idaho Batholith Fir-Spruce 6559 0.13
Fire_21_1988 1988 Middle Rockies Fir-Spruce 3113 0.17
Fire_22_1988 1988 Middle Rockies Other 29233 0.16
Fire_23_1988 1988 Middle Rockies Other 1363 0.41
Fire_24_1988 1988 Idaho Batholith Douglas-Fir 48136 0.31
Fire_25_1988 1988 Northern Rockies Douglas-Fir 8089 0.25
Fire_26_1988 1988 Middle Rockies Douglas-Fir 8588 0.56
Fire_27_1988 1988 Northern Rockies Douglas-Fir 11403 0.30
Fire_28_1988 1988 Northern Rockies Douglas-Fir 21854 0.25
Fire_29_1988 1988 Canadian Rockies Douglas-Fir 33844 0.27
Fire_30_1988 1988 Northern Rockies Douglas-Fir 1909 0.31
Fire_31_1988 1988 Middle Rockies Other 6282 0.47
Fire_32_1989 1989 Idaho Batholith Fir-Spruce 13334 0.15
Fire_33_1989 1989 Middle Rockies Lodegepole Pine 3298 0.37
Fire_34_1989 1989 Idaho Batholith Douglas-Fir 4928 0.38
Fire_35_1989 1989 Idaho Batholith Douglas-Fir 47680 0.02
Fire_36_1989 1989 Idaho Batholith Fir-Spruce 2486 0.30
Fire_37_1989 1989 Idaho Batholith Fir-Spruce 5566 0.30
Fire_38_1989 1989 Idaho Batholith Fir-Spruce 7443 0.28
Fire_39_1989 1989 Idaho Batholith Fir-Spruce 6786 0.26
Fire_40_1989 1989 Idaho Batholith Douglas-Fir 8733 0.12
Fire_41_1989 1989 Idaho Batholith Fir-Spruce 1615 0.49
Fire_42_1989 1989 Idaho Batholith Lodegepole Pine 2488 0.33
Fire_43_1989 1989 Idaho Batholith Fir-Spruce 3081 0.27
Fire_44_1990 1990 Middle Rockies Douglas-Fir 2535 0.45
Fire_45_1990 1990 Idaho Batholith Fir-Spruce 3418 0.53
Fire_46_1990 1990 Idaho Batholith Douglas-Fir 2249 0.49
Fire_47_1990 1990 Middle Rockies Lodegepole Pine 2763 0.41
Fire_48_1990 1990 Middle Rockies Other 13461 0.19
Fire_49_1991 1991 Middle Rockies Douglas-Fir 6978 0.31
Fire_50_1991 1991 Middle Rockies Douglas-Fir 3097 0.18
Fire_51_1991 1991 Middle Rockies Fir-Spruce 6995 0.14
Fire_52_1991 1991 Idaho Batholith Douglas-Fir 1186 0.55
Fire_53_1991 1991 Northern Rockies Fir-Spruce 7095 0.39
Fire_54_1991 1991 Northern Rockies Fir-Spruce 2478 0.51

2.4 Export Data

2.4.1 Final Cleanup for Export

# reformat and project
fires_export <- fires_select %>% 
  mutate(year = as.integer(year)) %>% 
  st_transform(., crs="EPSG:4326")

mtbs_export <- fires_mtbs %>% 
  mutate(year = as.integer(year)) %>% 
  st_transform(., crs="EPSG:4326")

2.4.2 Export

# st_write(fires_export, "data/fire_boundaries/", "fires_export.shp", driver = 'ESRI Shapefile')
# st_write(mtbs_export, "data/fire_boundaries/", "mbts_export.shp", driver = 'ESRI Shapefile')

2.5 Mapping Potential High-Severity Fires

2.5.1 Selected fires by year

# plot selected fire events by year
mapview(fires_select %>% mutate(year = as.factor(year)), zcol = "year",layer.name = "Fire Year") + mapview(st_union(eco_select),color = 'black', lwd = 1, alpha.regions = 0)

2.5.2 Selected fires by majority forest type

# plot selected fire events by majority forest type
mapview(fires_select %>% mutate(fire_foresttype = as.factor(fire_foresttype)), zcol = "fire_foresttype",layer.name = "Majority Forest Type") + mapview(st_union(eco_select),color = 'black', lwd = 1, alpha.regions = 0)

2.6 Mapping Final High-Severity Fires

2.6.1 Final Fires

# list of final fire events, after filtering to remove reburnt and managed areas
fire_names <- c("Fire_1_1988","Fire_2_1988","Fire_3_1988","Fire_4_1988","Fire_7_1988","Fire_9_1988","Fire_10_1988","Fire_11_1988","Fire_12_1988","Fire_13_1988","Fire_14_1988","Fire_15_1988","Fire_16_1988","Fire_18_1988","Fire_19_1988","Fire_20_1988","Fire_22_1988","Fire_23_1988","Fire_25_1988","Fire_26_1988","Fire_28_1988","Fire_29_1988","Fire_31_1988","Fire_32_1989","Fire_33_1989","Fire_35_1989","Fire_38_1989","Fire_41_1989","Fire_42_1989","Fire_48_1990","Fire_49_1991","Fire_50_1991","Fire_51_1991","Fire_54_1991")

2.6.2 Final Map

crs = 4388
states49 <- (us_states["NAME"])%>% st_transform(3488)
borders <- rnaturalearth::ne_states(c("united states of america", "canada")) %>% st_as_sf()
eco_border <- st_union(eco_select) %>% st_transform(crs = crs) 
conus_forestgroup <- rast('data/forest_type/conus_forestgroup.tif')

forestgroup_map <- crop(conus_forestgroup,st_buffer(eco_select,290000))%>% 
  as.factor() %>% 
  aggregate(fact=3,fun = "modal") %>% 
  project(.,eco_border)
## 
|---------|---------|---------|---------|
=========================================
                                          
forestgroup_df <- as.data.frame(forestgroup_map,xy = T) %>% 
  mutate(ftype = case_when(Red == 0 ~ "Unforested",
                           Red == 200 ~ "Douglas-Fir",
                           Red == 220 ~ "Ponderosa Pine",
                           Red == 280 ~ "Lodgepole Pine",
                           Red == 260 ~ "Fir-Spruce",
                           TRUE ~ "Other"))

forestgroup_df$ftype <- fct_relevel(forestgroup_df$ftype , c("Douglas-Fir","Lodgepole Pine","Fir-Spruce","Ponderosa Pine", "Other", "Unforested"))

palette <- brewer.pal(4,"YlGn")

inset <-ggplot() +
  geom_sf(data = states49 %>% st_transform(crs = crs), fill = NA,color=alpha("#525258"),lwd=.4) +
  geom_sf(data = st_union(eco_select) %>% st_transform(crs = crs),color = "black",fill="black",lwd=.4) +
  theme_bw() +
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())

mainplot <- ggplot() +
  geom_raster(data = forestgroup_df,aes(x,y,fill=as.factor(ftype)),alpha = .8)+
  scale_fill_manual(values = c("#C3E699","#68CE66","#238443","#8BB62F","#D8E895","white"))+ 
  geom_sf(data = borders %>% st_transform(crs = crs), fill = NA,color=alpha("black",0.14),lwd=.05) +
  geom_sf(data = fires_select %>% filter(Fire_ID %in% fire_names) %>% st_transform(crs = crs),fill = "black", color = "black") +
  geom_sf(data = st_union(eco_select) %>% st_transform(crs = crs), fill = NA, color=alpha("black",0.7),lwd=.4,linetype = "F1") +
  scale_x_continuous(limits = c(-119.4952-.5,-103.21640-1.68), expand = c(0, 0)) +
  scale_y_continuous(limits = c(41.89902-.5,49.00152+.5), expand = c(0, 0)) +
  labs(fill = "Forest Type", x = "Longitude", y = "Latitude") +
   theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())

ggdraw() +
  draw_plot(mainplot) +
  draw_plot(inset,
    height = 0.2,
    x = .39,
    y = 0.15
  )